%% 规划与图论Planning and graph theory
%% 1.整数规划Integer programming
% 下方给出整数规划的基本简单问题求解。
% (背包问题)总共有n件物品需要携带，每件物品价值为c(i),i=1,2,...,n，质量为
% a(i),i=1,2,...,n，要求最多携带b千克物品，则携带哪几件物品收益最大？
% 该问题为求解优化问题：max sum(c(i)*x(i)),i=1,2,...,n;s.t.sum(a(i)*x(i))<=b
% i=1,2,...,n;x(i)=0或1
clear,clc
n=10;% 10件物品
c=randi([1,10],1,10);% 价值
a=randi([1,10],1,10);% 质量
b=40;% 携带质量限制
f=@(x)sum(c.*x);% 目标函数
% 对于该问题，本节首先从连续角度思考解决方案，通过基础优化策略可以得知，将约
% 束条件通过罚函数法或增广拉格朗日函数法构造无约束优化问题，再利用无约束优化
% 算法可以求解问题，但整数规划中变量离散，构造罚函数需要较大的权值，容易使问
% 题产生不必要的变化，同时构造的函数有可能存在多个极小值情况，本节针对以上问
% 题，首先设计含有梯度下降的模拟退火求解方案。
% 对于连续型目标函数，可以设计如下程序：
clear,clc
lb=[-20,-20];
ub=[20,20];
n=length(lb);% 变量个数
% alpha0m=1/n*mean(ub-lb);% 初始步长
rand1=rand();
x0=rand1*ub+(1-rand1)*lb;% 迭代初始点
f=@(x)(x(1)-1)^2+(x(2)-2)^2;% 目标函数
time=100;% 进行100次求解
time1=30;% 每次求解进行30次
time2=1e3;% 求解过程迭代最大次数
result3=zeros(time,n+1);% 储存结果
Tm=100;% 退火初始温度即最大温度
dx=(1e-6)*(ub-lb);% 求梯度时维度微元设置
x=x0;
x1=x0;
% 迭代求解
for m=1:time
    x=lb+m/time*(ub-lb);% 迭代初始点
    result=zeros(time1,n+1);% 储存每次求解结果
    tic
    for k=1:time1
        T=100;% 退火温度初始化
        for n1=1:time2
            L1=f(x);
            k_un=-1+2*rand(1,n);% 随机方向
            z_un=k_un/norm(k_un);% 随机方向规范化
            c_u=zeros(1,n);% 该点梯度初始化
            for bianliang=1:n% 求梯度
                c_un1=x(bianliang)-dx(bianliang);
                c_un2=x(bianliang)+dx(bianliang);
                if x(bianliang)-dx(bianliang)<lb(bianliang)% 当要计算的梯度值小于变
                    % 量下限时处理
                    r_un=rand();
                    c_un1=r_un*lb(bianliang)+(1-r_un)*x(bianliang);
                elseif x(bianliang)+dx(bianliang)>ub(bianliang)% 当要计算的梯度值大
                    % 于变量下限时处理
                    r_un=rand();
                    c_un2=r_un*ub(bianliang)+(1-r_un)*x(bianliang);
                end
                xp=x;x2=x;
                xp(bianliang)=c_un1;x2(bianliang)=c_un2;
                c_u(bianliang)=(f(x2)-f(xp))/(2*dx(bianliang));
            end
            c_u=c_u/norm(c_u);% 确定该点最终梯度方向，当计算中未出现超过变量限制时为
            % 梯度方向，否则为次梯度方向。
            randchoice=rand();
            if(randchoice<0.2)
                z_un=-c_u;% 此处说明方向有40%概率为梯度方向，否则为随机方向
            end
            for i=1:n
                k1_un=x(i)+T/Tm*(ub(i)-lb(i))*z_un(i);
                if k1_un>=lb(i)&&k1_un<=ub(i)
                    x1(i)=k1_un;
                elseif k1_un<lb(i)
                    r_un=rand();
                    x1(i)=r_un*lb(i)+(1-r_un)*x(i);
                else
                    r_un=rand();
                    x1(i)=r_un*ub(i)+(1-r_un)*x(i);
                end
            end
            L2=f(x1);
            if L2<L1
                x=x1;
            else
                r_un=rand();
                p_un=exp(-abs(L2-L1)/T);
                if r_un<p_un
                    x=x1;
                end
            end
            T=(0.99*T);
        end
        if L2<L1
            result(k,:)=[L2,x1];
        else
            result(k,:)=[L1,x];
        end
    end
    result1=min(result(:,1));
    for m1=1:time1
        if result(m1,1)==result1
            result3(m,:)=result(m1,:);
        end
    end
    toc
end
% 利用以上程序即可多次求解优化问题，在概率上更有可能找到问题的解。
% 对于原始背包问题：
% 通过罚函数法构建拉格朗日函数
rho1=1e4;rho2=1e4;
L=@(x)-f(x)+rho1*sum(abs(x).*abs(x-1))+rho2*(max([sum(a.*x)-b,0]))^2;
lb=zeros(1,n)-1;% 求解下限
ub=1*ones(1,n)+1;% 求解上限
% 首先通过基本优化策略中线搜索方法设计求解函数，见文件linesearch_optim.m
[f_min,x]=linesearch_optim(L,lb,ub);
% 之后通过含有梯度下降的模拟退火设计求解函数，见文件Grad_Simulated_annealing.m
result=Grad_Simulated_annealing(L,lb,ub);
% 由分析结果可知，连续型优化策略可以找到较好的点，但求解结果不稳定，比较容易
% 收敛到局部最小解。
% matlab基于问题的求解器可以解决小型整数型规划问题，则对于背包问题，有：
x=optimvar('x',n,'Type','integer','LowerBound',0,'UpperBound',1);% 声明变量
% 类型
prob=optimproblem("ObjectiveSense","min");
prob.Objective=-c*x;% 目标函数
prob.Constraints.con=a*x-b<=0;
[sol,fval,flag]=solve(prob);
% 可见针对离散型数据问题需要用离散型解决方案进行解决，连续型的解决策略并不适
% 合求解该类问题。本节讨论几种0-1整数规划(背包问题)的解决方案。
% 枚举法解决方案：整数规划相比较连续型，其只有有限个解，因此可以设计具有过滤
% 条件的枚举法解决整数规划问题。
% (该问题程序可见文件：enumeration0_1.m)
clear,clc
n=10;% 变量个数
c=[1,6,9,8,3,4,7,5,7,1];
b=40;
a=[2,4,10,2,6,5,4,6,9,6];
f=@(x)-sum(c.*x);% 目标函数
con=@(x)sum(a.*x)-b<=0;% 约束条件
x=zeros(1,n);% 初始计算点
x1=zeros(1,n);
result0=100;% 初始点函数值，给一个较大值，以便迭代正常更新
nm=floor(30*log(10)/log(2));% 自变量个数阈值
if(n<nm)
    p=2^n;
    for i=0:p-1
        xstr=dec2bin(i);
        xstr1=mat2cell(xstr,1,ones(1,numel(xstr)));
        xstr2=strjoin(xstr1,' ');
        xstr2=[' ',xstr2];
        len=numel(xstr);
        if(len~=n)
            xstr0=zeros(1,n-len);
            xstr2=[num2str(xstr0),xstr2];
        end
        x=str2num(xstr2);% 进行数值转换操作
        if(con(x)==1)
            result1=f(x);
            if(result1<result0(1))
                result0=result1;
                x1=x;
            elseif(result1==result0(1))
                result0=[result0;result1];
                x1=[x1;x];
            end
        end
    end
end
% 可见当变量个数不大时，枚举法具有较高的运行速度，且必定可以找到0-1规划问题
% 的最优解。
% 当约束条件为整数时，可以放松条件先求解连续问题的最优解，再在连续问题最优解
% 附近寻找合适的取值。可以采取分支定界，割平面法进行求解。
% 当取值范围为0,1,2;0,1,2,3等时，可以利用除基取余的方法进行枚举计算，如下所示。
ne=35469;basic=3;
test=[];ne1=ne;
while ne1>=basic
    rem0=mod(ne1,basic);
    ne1=round((ne1-rem0)/basic);
    test=[test,rem0];
end
test=[test,ne1];test1=fliplr(test);
%% 2.典型非线性规划Typical nonlinear programming
% 本节通过一飞行管理问题建模分析不同非线性优化算法。
% 在约10000m高空某160km*160km正方形区域内，经常有若干架飞机水平飞行。区域内
% 每架飞机的位置和速度向量由计算机记录。当一架飞机进入该区域到达区域边缘时，
% 记录其数据后，要立即计算并判断是否会余区域内飞机发生碰撞，若会发生碰撞，则
% 应计算如何调整各架飞机飞行的方向角，以避免碰撞。现有如下条件需要满足：
% (1)不碰撞标准为任意两架飞机距离大于8km；
% (2)飞机飞行角调整幅度不超过30度；
% (3)所有飞机飞行速度均为800km/h；
% (4)进入该区域飞机在到达区域边缘时，与区域内飞机距离应在60km以上；
% (5)最多需考虑6架飞机；
% (6)不必考虑飞机离开此区域后状况。
% 建模计算下方数据模型，要求每架飞机飞行方向角调整幅度尽可能小。
% 设该区域4个顶点为：(0,0),(160,0),(160,160),(0,160)。
% 飞机编号，横坐标，纵坐标，方向角见flight.mat
clear,clc
load flight.mat % 加载数据
n=6;% 飞机架数
theta=zeros(1,n);% 每架飞机调整角初始化
dt=1e-4;% 时间微元
% 本节首先考虑对约束条件不碰撞和区域边缘限制进行软化，通过设计处罚，以使优化
% 结果满足条件限制。
tm=sqrt(2*160^2)/800;% 最大运行时间
v=800;% 速度
flight1=flight;
flight1(:,4)=flight1(:,4)/360*2*pi;% 将角度转换为弧度制
d_crash=zeros(1,n);% 区域内碰撞处罚参数
d_crashm=zeros(2,n);% 飞机到达区域边缘时距离限制处罚参数
for i=0:dt:tm
    for i1=1:n % 模拟飞行过程
        if(flight1(i1,2)>=0&&flight1(i1,2)<=160&&flight1(i1,3)>=0&&flight1(i1,3)<=160)
            flight1(i1,2)=flight1(i1,2)+v*dt*cos(flight1(i1,4)+theta(i1));
            flight1(i1,3)=flight1(i1,3)+v*dt*sin(flight1(i1,4)+theta(i1));
        end
    end
    for i1=1:n % 计算区域内飞机间距离
        if(flight1(i1,2)>0&&flight1(i1,2)<160&&flight1(i1,3)>0&&flight1(i1,3)<160)
            d1=0;
            for i2=1:n
                if(i1~=i2)
                    if(flight1(i2,2)>0&&flight1(i2,2)<160&&flight1(i2,3)>0&& ...
                            flight1(i2,3)<160)
                        d0=sqrt((flight1(i1,2)-flight1(i2,2))^2+(flight1(i1,3)- ...
                            flight1(i2,3))^2);
                        if(d0<8)
                            d1=d1+1/(d0+eps)*dt/tm;
                        end
                    end
                end
            end
            d_crash(i1)=d_crash(i1)+d1;
        end
    end
    for i1=1:n % 计算飞机在区域边缘时情况
        if(flight1(i1,2)<=0||flight1(i1,2)>=160||flight1(i1,3)<=0||flight1(i1,3)>=160)
            if(d_crashm(1,i1)==0)
                d1=0;
                for i2=1:n
                    if(i1~=i2)
                        if(flight1(i2,2)>0&&flight1(i2,2)<160&&flight1(i2,3)>0&& ...
                                flight1(i2,3)<160)
                            d0=sqrt((flight1(i1,2)-flight1(i2,2))^2+(flight1(i1,3)- ...
                                flight1(i2,3))^2);
                            if(d0<60)
                                d1=d1+60/(d0+eps);
                            end
                        end
                    end
                end
                d_crashm(2,i1)=d_crashm(2,i1)+d1;
                d_crashm(1,i1)=1;
            end
        end
    end
end
% 通过加入处罚因子即可设计目标函数：
rho1=100;rho2=100;% 处罚因子
f_min=sum(abs(theta))/pi*180+rho1*sum(d_crash)+rho2*sum(d_crashm(2,:));% 目标函数设计
% 设计函数见文件:example_flight.m
% 接下来考虑利用不同优化方法寻找最优解。
f=@(theta)example_flight(theta,flight);
lb=-pi/6*ones(1,n);
ub=pi/6*ones(1,n);% 范围限制
% 梯度下降线搜索法：linesearch_optim.m
[f_min,theta]=linesearch_optim(f,lb,ub);
theta1=theta/pi*180;
% 利用线搜索方法找到较好解有：4.78 	-1.76 	-3.98 	-8.60 	3.46 	-0.47 
% 7.73 	-0.06 	-7.08 	-7.68 	5.12 	0.05 
% 可见利用线搜索梯度下降结果具有不稳定性。
% 含有梯度下降的模拟退火法：Grad_Simulated_annealing.m
result=Grad_Simulated_annealing(f,lb,ub);
theta1=result(:,2:end)/pi*180;
% 较好解有：7.2682 	-0.0036 	-4.7274 	-4.4281 	4.2359 	-0.8991 
% 6.7162 	0.0023 	-3.1478 	-2.7880 	2.6395 	-0.0044 
% 8.8923 	4.6709 	-5.3984 	-4.6604 	7.1726 	-0.1620 
% 由于牛顿类算法和信赖域算法需要用到目标函数二阶导数信息，过程较为复杂，不适
% 合本题求解。
%% 3.图论与网络模型Graph theory and network model
%%  3.1.最短路算法
% 设图G是赋权图，T为图G中的一条路，则称T的各边权之和为路T的长度。对于连通图
% (任意两顶点之间存在道路)G的两个顶点之间的路可能不止一条，其中权值之和最小的
% 路称为两顶点间最短路。
% (有向图求最短路问题)已知某人要从v1出发前往v8，求最短路。
% 该图邻接矩阵如下，其中，逆行和未连接点间权值为1e4：
clear,clc
w=[0,6,3,1,1e4,1e4,1e4,1e4,1e4
    1e4,0,1e4,1e4,1,1e4,1e4,1e4,1e4
    1e4,2,0,2,1e4,1e4,1e4,1e4,1e4
    1e4,1e4,1e4,0,1e4,10,1e4,1e4,1e4
    1e4,1e4,1e4,6,0,4,3,6,1e4
    1e4,1e4,1e4,1e4,10,0,2,1e4,1e4
    1e4,1e4,1e4,1e4,1e4,1e4,0,4,1e4
    1e4,1e4,1e4,1e4,1e4,1e4,1e4,0,1e4
    1e4,1e4,1e4,1e4,2,1e4,1e4,3,0];
% 可绘制出该有向图：
E=[1,2,6;1,3,3;1,4,1;2,5,1;3,2,2;3,4,2;4,6,10;5,4,6
    5,6,4;5,7,3;5,8,6;6,5,10;6,7,2;7,8,4;9,5,2;9,8,3];
G=digraph(E(:,1),E(:,2),E(:,3));
% 求最短路可看作一0-1规划问题，目标函数为：min sum(sum(w.*x))，当端点i到j在
% 最短路上时，令x(i,j)=1。
% 则有如下目标函数与限制条件：
f=@(x)sum(sum(w.*x));
con1=@(x)sum(x(1,:))==1;% 从起始点1出发
con2=@(x)sum(x(:,8))==1;% 最终到达点8
con3=@(x)sum(x(:,1))==0;% 从1出发一定到达其他点
con4=@(x)prod([sum(x(:,2:7)),sum(x(:,9))]==[sum(x(2:7,:),2);sum(x(9,:))]');% 从某点
% 出发必定有其他点到达该点
% 利用传统0-1枚举法求解问题需要2^81次运算，显然不可能完成。仔细思考问题的约
% 束条件，x每一横行代表从某点出发前往令一个点，因此x为1的点只有出现在w对应行
% 不为0或1e4的位置，因此计算次数最多为2^16次，计算机可以快速求解出问题。
n=16;% 共进行2^16次运算
x0=zeros(9,9);% 参考位置点
for i=1:9
    for i1=1:9
        if(w(i,i1)~=0&&w(i,i1)<1e4)
            x0(i,i1)=1;
        end
    end
end
x1=[];
result0=100;% 初始点函数值，给一个较大值，以便迭代正常更新
nm=floor(30*log(10)/log(2));% 自变量个数阈值
if(n<nm)
    p=2^n;
    for i=1:p
        % 利用除基取余转换为2进制
        ne=i;basic=2;
        test=[];ne1=ne;
        while ne1>=basic
            rem0=mod(ne1,basic);
            ne1=round((ne1-rem0)/basic);
            test=[test,rem0];
        end
        test=[test,ne1];test1=fliplr(test);
        len=length(test1);
        if(len<=n-1)
            xp0=zeros(1,n-len);
            xp=[xp0,test1];
        else
            xp=test1;
        end
        p0=0;x=zeros(9,9);
        for i1=1:9
            for i2=1:9
                if(x0(i1,i2)==1)
                    p0=p0+1;
                    x(i1,i2)=xp(p0);
                end
            end
        end
        k1=sum(x(:,2:7));k2=sum(x(:,9));
        k3=sum(x(2:7,:),2);k4=sum(x(9,:));
        p0=[k1,k2]==[k3;k4]';
        p1=prod(p0);
        if(con1(x)&&con2(x)&&con3(x)&&p1)
            result1=f(x);
            if(result1<result0(1))
                result0=result1;
                x1=x;
            elseif(result1==result0(1))
                result0=[result0;result1];
                x1=[x1;x];
            end
        end
    end
end
% 求解的最短路长度为12，为：1-3-2-5-8
% 当图被绘制出来时，matlab含有内置的求解最短路程序:
[path,d]=shortestpath(G,1,8,'Method','positive');% 其中path为路径，d为长度
% 由以上求解过程可知，求最短路时可行解的范围是有限的，甚至远远小于上方0-1规
% 划使用到的求解次数。本节介绍求解最短路的Dijkstra算法，Floyd算法。
% Dijkstra算法：
% Dijkstra算法依据是：最短路是一条路，最短路上任意一子段也是最短路。
% 按距固定起点u0从近到远为顺序，依次求得u0到图G各顶点的最短路和距离，直至某
% 个顶点v0或图G所有顶点。
clear,clc
w=[0,6,3,1,1e4,1e4,1e4,1e4,1e4
    1e4,0,1e4,1e4,1,1e4,1e4,1e4,1e4
    1e4,2,0,2,1e4,1e4,1e4,1e4,1e4
    1e4,1e4,1e4,0,1e4,10,1e4,1e4,1e4
    1e4,1e4,1e4,6,0,4,3,6,1e4
    1e4,1e4,1e4,1e4,10,0,2,1e4,1e4
    1e4,1e4,1e4,1e4,1e4,1e4,0,4,1e4
    1e4,1e4,1e4,1e4,1e4,1e4,1e4,0,1e4
    1e4,1e4,1e4,1e4,2,1e4,1e4,3,0];
% 可绘制出该有向图：
E=[1,2,6;1,3,3;1,4,1;2,5,1;3,2,2;3,4,2;4,6,10;5,4,6
    5,6,4;5,7,3;5,8,6;6,5,10;6,7,2;7,8,4;9,5,2;9,8,3];
G=digraph(E(:,1),E(:,2),E(:,3));
numsize=size(w,2);% 节点个数
start=1;% 起始节点
s=start;% 具有永久标号的顶点集初始化
dist=w(start,:);% 距离矩阵初始化
p=zeros(1,numsize);% 前驱矩阵初始化
vs=zeros(1,numsize);% 未激活顶点初始化
for i=1:numsize
    if(dist(i)==0||dist(i)>=1e4)
        p(i)=-1;
    else
        p(i)=1;
    end
    if(i~=start)
        vs(i)=i;
    end
end
dist1=zeros(1,numsize);
% 迭代，通过找最小，加入s战队，借东风三步逐步寻找最短路
while sum(vs)~=0 % 当还存在未激活点时
    % 找最小
    t=max(dist1);d0=[];
    for i=1:numsize
        if(dist(i)>t)
            d0=[d0,dist(i)];
        end
    end
    d1=min(d0);
    for i=1:numsize
        if(dist(i)==d1)
            dist1(i)=dist(i);
            t0=i; % 最小节点
        end
    end
    % 加入s战队
    s=[s,t0];vs(t0)=0;
    % 借东风
    dist2=w(t0,:);d2=[];
    for i=1:numsize
        if(dist2(i)~=0&&dist2(i)<1e4)
            d2=[d2,[i;dist2(i)]];
        end
    end
    if(size(d2,2)~=0)
        for i=1:length(d2(1,:))
            if(dist(t0)+d2(2,i)<dist(d2(1,i)))
                dist(d2(1,i))=dist(t0)+d2(2,i);
                p(d2(1,i))=t0;
            end
        end
    end
end
% 利用前驱矩阵p即可找到最短路，最短路长度在dist或dist1中。
% Dijkstra算法要求所有边或弧的权都是非负的。而Floyd算法允许赋权图中包含负权
% 的边或弧。
% Floyd算法：
% 求距离矩阵，路径矩阵：
b=w;% 获得邻接矩阵
r=zeros(numsize);
for i=1:numsize
    for i1=1:numsize
        for i2=1:numsize
            if(b(i1,i)+b(i,i2)<b(i1,i2))
                b(i1,i2)=b(i1,i)+b(i,i2);
                r(i1,i2)=i;
            end
        end
    end
end
% 获得矩阵b元素即为某两点间最短距离。
% 矩阵r中，若r(i,j)=p，则点vp是顶点vi到vj最短路的中间点。
%% 3.2.最小生成树
% 在赋权图G中，边权之和最小的生成树称为G的最小生成树。
% Kruskal算法：
% 每次将一条权最小的边加入子图T中，并保证不形成圈。
clear,clc
w=[0,20,33,63,15,30;20,0,20,50,25,40;33,20,0,30,18,33;63,50,30,0,48,63
    15,25,18,48,0,15;30,40,33,63,15,0];% 邻接矩阵
s=[];% 最小生成数顶点集
s1=[];s2=[];
t0=0;t2=[];
numsize=size(w,2);
p=cumsum(ones(1,numsize));
while sum(p~=0)
    t1=min(w(w>t0));
    [i,j]=find(w==t1);
    for i1=1:length(j)
        p5=[i(i1),j(i1)];
        p6=0;
        if(size(s2,2)~=0)
            for i2=1:length(s2(:,1))
                p6=p6+prod(s2(i2,:)==p5);
            end
        end
        if(p6==0)
            s1=[s1;p5];
            t2=[t2;t1];
        end
    end
    t0=t1;
    for i=1:numsize
        for i1=1:length(s1(:,1))
            if(p(i)==s1(i1,1))
                p(i)=0;
            end
        end
    end
    p1=zeros(1,numsize);
    for i=1:numsize
        p2=0;
        for i1=1:length(s1(:,1))
            if(s1(i1,1)==i)
                p2=p2+1;
            end
        end
        p1(i)=p2;
    end
    s2=[];
    for i=1:numsize
        s3=[];
        if(p1(i)>=2)
            p4=0;p3=0;
            for i1=1:length(s1(:,1))
                if(s1(i1,1)==i)
                    p4=p3;
                    p3=s1(i1,2);
                    if(p3~=p4)
                        s3=[s3,p3];
                    end
                end
            end
        end
        if(size(s3,2)~=0)
            s3=nchoosek(s3,2);
            s2=[s2;s3];
        end
    end
    s2=[s2;fliplr(s2)];
end
s=[s1,t2];
% 利用如上程序算法可以找到最小生成树所在的边和点，算法可以保证不会出现3点构
% 成圈，但可能出现3点以上构成的圈，此时只需将圈中较大的边断开即可。
% matlab中含有内置函数可以求最小生成树：
sv=cellstr(strcat('v',int2str((1:numsize)')));% 建立顶点元胞
G=graph(w,sv,'upper');% 绘制图
T=minspantree(G,'Method','sparse');% Kruskal算法
pl=plot(G,'EdgeLabel',G.Edges.Weight);% 画图
highlight(pl,T)
%% 3.3.着色问题、最大流和旅行商问题
% 着色问题：对图G所有顶点着色，要求相邻顶点颜色不同，求使用到最小的颜色数。
% 该问题可转化为0-1规划。
clear,clc
w=[0,1,1e4,1,1;1,0,1,1,1e4;1e4,1,0,1,1;1,1,1,0,1e4;1,1e4,1,1e4,0];% 邻接矩阵
% 当两顶点间有连接时为1，否则为1e4.
numsize=size(w,2);
w1=zeros(numsize);
for i=1:numsize
    for i1=1:numsize
        if(w(i,i1)==1)
            w1(i,i1)=1;
        end
    end
end
[ni,nj]=find(w1);% 边的顶点编号
deg=sum(w1);Km=max(deg);% 计算顶点的度及最大度
prob=optimproblem;% 建立0-1优化问题模型求解着色问题
x=optimvar('x',numsize,Km+1,'Type','integer','LowerBound',0,'UpperBound',1);% 建立
% 0-1规划变量
y=optimvar('y');prob.Objective=y;
prob.Constraints.con1=sum(x,2)==1;
prob.Constraints.con2=x(ni,:)+x(nj,:)<=1;
prob.Constraints.con3=x*(1:Km+1)'<=y;
[sol,fval,flag,out]=solve(prob);
% 最大流与最小费用流问题：求解该类问题可使用Ford-Fulkerson等方法，限于文章篇
% 幅，本文不再过多介绍。最大流和最小费用流问题可以写为线性规划形式，matlab中
% 含有求最大流的函数。
clear,clc
NN=cellstr(strcat('v',int2str((2:5)')));% 发点为1，终点为6，构造中间节点
NN=[{'vs'},NN(:)',{'vt'}];
L={'vs','v2',5,3;'vs','v3',3,6;'v2','v4',2,8;'v3','v2',1,2;'v3','v5',4,2
    'v4','v3',1,1;'v4','v5',3,4;'v4','vt',2,10;'v5','vt',5,2};% 构造网络
% 其中第一个数字为两点间容量，第二个数字为两点间单位运输费用。
G=digraph;G=addnode(G,NN);
G1=addedge(G,L(:,1),L(:,2),cell2mat(L(:,3)));% 导出流量到图中
[M,F]=maxflow(G1,'vs','vt');% 求最大流
G2=addedge(G,L(:,1),L(:,2),cell2mat(L(:,4)));% 导出单位费用到图中
c=full(adjacency(G1,"weighted"));% 导出容量矩阵
b=full(adjacency(G2,"weighted"));% 导出费用矩阵
prob=optimproblem;
f=optimvar('f',6,6,'LowerBound',0);
prob.Objective=sum(sum(b.*f));
con1=[sum(f(1,:))==M;sum(f(:,2:end-1))'==sum(f(2:end-1,:),2)
    sum(f(:,end))==M];
prob.Constraints.con1=con1;
prob.Constraints.con2=f<=c;
[sol,fval,flag,out]=solve(prob);
% 旅行商问题：一名推销员准备前往若干城市推销产品，然后回到他的出发地。如何为
% 他设计最短的旅行路线(从驻地出发，经过每个城市恰好一次，最后返回驻地)。可以
% 通过改良圈算法求得较好解，也可以通过0-1规划求解。
clear,clc
w=[0,8,5,9,12,14,12,16,17,22;8,0,9,15,17,8,11,18,14,22
    5,9,0,7,9,11,7,12,12,17;9,15,7,0,3,17,10,7,15,18
    12,17,9,3,0,8,10,6,15,15;14,8,11,17,8,0,9,14,8,16
    12,11,7,10,10,9,0,8,6,11;16,18,12,7,6,14,8,0,11,11
    17,14,12,15,15,8,6,11,0,10;22,22,17,18,15,16,11,11,10,0];% 邻接矩阵
numsize=size(w,2);
w1=w;
w1(1:numsize+1:end)=1e4;% 将对角线元素换为充分大
prob=optimproblem;
x=optimvar('x',numsize,numsize,'Type','integer','LowerBound',0,'UpperBound',1);
u=optimvar('u',numsize,'LowerBound',0);% 序号变量
prob.Objective=sum(sum(w1.*x));
prob.Constraints.con1=[sum(x,2)==1;sum(x,1)'==1;u(1)==0];
con2=[1<=u(2:end);u(2:end)<=14];
for i=1:numsize
    for j=2:numsize
        con2=[con2;u(i)-u(j)+numsize*x(i,j)<=numsize-1];
    end
end
prob.Constraints.con2=con2;
[sol,fval,flag]=solve(prob);
%% 4.多目标规划与目标规划Multi-objective programming and objective programming
% 多目标规划相比较常规规划问题，目标函数具有多个，一般可通过以下几种常用方法
% 进行处理：
% 首先需要对每一个目标函数进行无量纲化处理，数量级归一化处理。
% (1)对目标函数线性加权求和作为单目标函数；
% (2)选取一项主要的目标函数作为求解目标，对其他目标函数各取参数e，将其他目标
% 函数作为约束条件；
% (3)求解每一个目标函数下对应最优解，构造每个目标与理想值差的加权平方和作为作
% 为评价函数；
% (4)确定每一个目标函数的优先级，先求最高级目标函数，将最高级目标函数等于最优
% 作为约束，求第二级目标函数，以此类推。
% 对应决策类优化问题，其一般对目标函数有明确的目标值，对于具体的决策问题，一
% 般情况下要明确以下概念并建立目标规划模型：
% (1) 正负偏差变量：设fi为第i个目标函数，它的正偏差变量di(+)=max{fi-di0,0}
% 负偏差变量di(-)=-min{fi-di0,0}，其中di0为fi的目标值；
% (2) 绝对刚性约束：只优化问题中必须严格满足的约束条件；
% (3) 目标约束：目标约束是目标规划特有的，其约束右端项是要达到的目标，左端为对
% 应目标函数加入正、负偏差；
% (4) 优先因子与权系数：要求第一位达到的目标赋予优先因子P1，次位目标赋予优先因
% 子P2，依次类推；
% (5) 目标规划的目标函数一般有以下三类：要求恰好达到目标值，即正负偏差变量尽
% 可能小；要求不超过目标值，即正偏差变量尽可能小；要求超过目标值，即负偏差变
% 量尽可能小。